####Model results####
##a) fiscal: revenue (table X32, models A111-A114)
r3b_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c("sa_revenue","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1))
r3b_m1g_t1_sc_cse <- data.frame(cluster.se(r3b_m1g_t1_sc, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
r3b_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c("sa_revenue","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0))
r3b_m1g_t1_sf_cse <- data.frame(cluster.se(r3b_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
r3b_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_revenue","included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r3b_m1d_t1_sb_cse <- data.frame(cluster.se(r3b_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
r3b_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_revenue*included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r3b_m2d_t1_sb_cse <- data.frame(cluster.se(r3b_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(r3b_m1g_t1_sc, r3b_m1g_t1_sf, r3b_m1d_t1_sb, r3b_m2d_t1_sb, se=c(r3b_m1g_t1_sc_cse, r3b_m1g_t1_sf_cse, r3b_m1d_t1_sb_cse, r3b_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A111)", "Min. (model A112)","Maj./min. dyad (model A113)","Maj./min. dyad (model A114)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X32. Additional results: alternative measurement of territorial autonomy - fiscal revenue.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Included","Included/excluded","Revenue x included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Revenue","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex32.txt")
stargazer(r3b_m1g_t1_sc, r3b_m1g_t1_sf, r3b_m1d_t1_sb, r3b_m2d_t1_sb, se=c(r3b_m1g_t1_sc_cse, r3b_m1g_t1_sf_cse, r3b_m1d_t1_sb_cse, r3b_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A111)", "Min. (model A112)","Maj./min. dyad (model A113)","Maj./min. dyad (model A114)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X32. Additional results: alternative measurement of territorial autonomy - fiscal revenue.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Included","Included/excluded","Revenue x included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Revenue","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex32.html")
##b) fiscal: expenditure (table X33, models A115-A118)
r3c_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c("sa_expenditure","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1))
r3c_m1g_t1_sc_cse <- data.frame(cluster.se(r3c_m1g_t1_sc, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
r3c_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c("sa_expenditure","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0))
r3c_m1g_t1_sf_cse <- data.frame(cluster.se(r3c_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
r3c_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_expenditure","included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r3c_m1d_t1_sb_cse <- data.frame(cluster.se(r3c_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
r3c_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("sa_expenditure*included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r3c_m2d_t1_sb_cse <- data.frame(cluster.se(r3c_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(r3c_m1g_t1_sc, r3c_m1g_t1_sf, r3c_m1d_t1_sb, r3c_m2d_t1_sb, se=c(r3c_m1g_t1_sc_cse, r3c_m1g_t1_sf_cse, r3c_m1d_t1_sb_cse, r3c_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A115)", "Min. (model A116)","Maj./min. dyad (model A117)","Maj./min. dyad (model A118)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X33. Additional results: alternative measurement of territorial autonomy - fiscal expenditure.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Included","Included/excluded","Expenditure x included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Expenditure","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex33.txt")
stargazer(r3c_m1g_t1_sc, r3c_m1g_t1_sf, r3c_m1d_t1_sb, r3c_m2d_t1_sb, se=c(r3c_m1g_t1_sc_cse, r3c_m1g_t1_sf_cse, r3c_m1d_t1_sb_cse, r3c_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A115)", "Min. (model A116)","Maj./min. dyad (model A117)","Maj./min. dyad (model A118)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X33. Additional results: alternative measurement of territorial autonomy - fiscal expenditure.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Included","Included/excluded","Expenditure x included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Expenditure","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex33.html")
##c) RAI self-rule (table X34, models A119-A122)
r3d_m1g_t1_sc <- glm(as.formula(paste("cw_event_g", paste(c("rai_autonomy","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1))
r3d_m1g_t1_sc_cse <- data.frame(cluster.se(r3d_m1g_t1_sc, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode))[,2])
r3d_m1g_t1_sf <- glm(as.formula(paste("cw_event_g", paste(c("rai_autonomy","included", group_nat_vars, group_unit_vars, unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cw_event_g_peaceyears_l1","I(cw_event_g_peaceyears_l1^2)","I(cw_event_g_peaceyears_l1^3)","cw_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0))
r3d_m1g_t1_sf_cse <- data.frame(cluster.se(r3d_m1g_t1_sf, as.factor(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode))[,2])
r3d_m1d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("rai_autonomy","included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r3d_m1d_t1_sb_cse <- data.frame(cluster.se(r3d_m1d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
r3d_m2d_t1_sb <- glm(as.formula(paste("cv_event", paste(c("rai_autonomy*included_excluded","excluded_excluded", dyad_nat_vars, dyad_unit_vars, unit_vars, nat_vars,"as.factor(region)","as.factor(year)","cv_event_peaceyears_l1","I(cv_event_peaceyears_l1^2)","I(cv_event_peaceyears_l1^3)","cv_event_any_slag_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1))
r3d_m2d_t1_sb_cse <- data.frame(cluster.se(r3d_m2d_t1_sb, as.factor(subset(main_dyad, subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode))[, 2])
stargazer(r3d_m1g_t1_sc, r3d_m1g_t1_sf, r3d_m1d_t1_sb, r3d_m2d_t1_sb, se=c(r3d_m1g_t1_sc_cse, r3d_m1g_t1_sf_cse, r3d_m1d_t1_sb_cse, r3d_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A119)", "Min. (model A120)","Maj./min. dyad (model A121)","Maj./min. dyad (model A122)"), omit=c("cowcode|region|factor"), type="text", order=vars.order_gd, title="Table X34. Additional results: alternative measurement of territorial autonomy - RAI self-rule index.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Included","Included/excluded","RAI self-rule x included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","RAI self-rule","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex34.txt")
stargazer(r3d_m1g_t1_sc, r3d_m1g_t1_sf, r3d_m1d_t1_sb, r3d_m2d_t1_sb, se=c(r3d_m1g_t1_sc_cse, r3d_m1g_t1_sf_cse, r3d_m1d_t1_sb_cse, r3d_m2d_t1_sb_cse), dep.var.labels.include = T, dep.var.labels = c("Civil violence (group)","Communal violence","Communal violence"), model.numbers =F, column.labels = c("Maj. (model A119)", "Min. (model A120)","Maj./min. dyad (model A121)","Maj./min. dyad (model A122)"), omit=c("cowcode|region|factor"), type="html", order=vars.order_gd, title="Table X34. Additional results: alternative measurement of territorial autonomy - RAI self-rule index.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses; maj. = second-order majority; min. = second-order minority.","The dependent variable is a binary variable equal to one if there is at least one instance of civil/communal violence involving a group/dyad in a given unit.","Region- and year-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Included","Included/excluded","RAI self-rule x included/excluded","Excluded/excluded","Relative size (state)","Relative size (state, mean)","Relative size (state, diff.)","Relative size (unit)","Relative size (unit, mean)","Relative size (unit, diff.)","Asymmetry", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","RAI self-rule","Civil violence peace years","Civil violence peace years 2","Civil violence peace years 3","Civil violence spatial lag", "Communal violence peace years", "Communal violence peace years 2", "Communal violence peace years 3","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex34.html")

####Figure A18####
##a) fiscal (revenue) (civil)
me_r3b_m1g_t1_sc <- margins_summary(r3b_m1g_t1_sc, variables = "sa_revenue", vcov = cluster.vcov(r3b_m1g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode)))
me_r3b_m1g_t1_sc$term <- "second-order maj."
me_r3b_m1g_t1_sf <- margins_summary(r3b_m1g_t1_sf, variables = "sa_revenue", vcov = cluster.vcov(r3b_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_r3b_m1g_t1_sf$term <- "second-order min."
me_r3b_m1g <- rbind.fill(me_r3b_m1g_t1_sc, me_r3b_m1g_t1_sf)
me_r3b_m1g$term <- as.factor(me_r3b_m1g$term)
me_r3b_m1g$term = factor(me_r3b_m1g$term,levels(me_r3b_m1g$term)[c(2,1)])
me_r3b_m1g$model <- "fiscal: revenue"
me_r3b_m1g$estimate <- me_r3b_m1g$AME
me_r3b_m1g$conf.low <- me_r3b_m1g$lower
me_r3b_m1g$conf.high <- me_r3b_m1g$upper
##b) fiscal (revenue) (communal)
me_r3b_m1d_t1_sb <- margins_summary(r3b_m1d_t1_sb, variables = "sa_revenue", vcov = cluster.vcov(r3b_m1d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r3b_m1d_t1_sb$term <- "all"
me_r3b_m1d_t1_sb$model <- "fiscal: revenue"
me_r3b_m2d_t1_sb <- margins_summary(r3b_m2d_t1_sb, variables = "sa_revenue", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(r3b_m2d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r3b_m2d_t1_sb$term <- ifelse(me_r3b_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_r3b_m2d_t1_sb$model <- "fiscal: revenue"
me_r3b_m13d <- rbind.fill(me_r3b_m1d_t1_sb, me_r3b_m2d_t1_sb)
me_r3b_m13d <- subset(me_r3b_m13d, term != "other")
me_r3b_m13d$term <- as.factor(me_r3b_m13d$term)
me_r3b_m13d$estimate <- me_r3b_m13d$AME
me_r3b_m13d$conf.low <- me_r3b_m13d$lower
me_r3b_m13d$conf.high <- me_r3b_m13d$upper
##c) fiscal (expenditure) (civil)
me_r3c_m1g_t1_sc <- margins_summary(r3c_m1g_t1_sc, variables = "sa_expenditure", vcov = cluster.vcov(r3c_m1g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode)))
me_r3c_m1g_t1_sc$term <- "second-order maj."
me_r3c_m1g_t1_sf <- margins_summary(r3c_m1g_t1_sf, variables = "sa_expenditure", vcov = cluster.vcov(r3c_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_r3c_m1g_t1_sf$term <- "second-order min."
me_r3c_m1g <- rbind.fill(me_r3c_m1g_t1_sc, me_r3c_m1g_t1_sf)
me_r3c_m1g$term <- as.factor(me_r3c_m1g$term)
me_r3c_m1g$term = factor(me_r3c_m1g$term,levels(me_r3c_m1g$term)[c(2,1)])
me_r3c_m1g$model <- "fiscal: expenditure"
me_r3c_m1g$estimate <- me_r3c_m1g$AME
me_r3c_m1g$conf.low <- me_r3c_m1g$lower
me_r3c_m1g$conf.high <- me_r3c_m1g$upper
##d) fiscal (expenditure) (communal)
me_r3c_m1d_t1_sb <- margins_summary(r3c_m1d_t1_sb, variables = "sa_expenditure", vcov = cluster.vcov(r3c_m1d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r3c_m1d_t1_sb$term <- "all"
me_r3c_m1d_t1_sb$model <- "fiscal: expenditure"
me_r3c_m2d_t1_sb <- margins_summary(r3c_m2d_t1_sb, variables = "sa_expenditure", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(r3c_m2d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r3c_m2d_t1_sb$term <- ifelse(me_r3c_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_r3c_m2d_t1_sb$model <- "fiscal: expenditure"
me_r3c_m13d <- rbind.fill(me_r3c_m1d_t1_sb, me_r3c_m2d_t1_sb)
me_r3c_m13d <- subset(me_r3c_m13d, term != "other")
me_r3c_m13d$term <- as.factor(me_r3c_m13d$term)
me_r3c_m13d$estimate <- me_r3c_m13d$AME
me_r3c_m13d$conf.low <- me_r3c_m13d$lower
me_r3c_m13d$conf.high <- me_r3c_m13d$upper
##e) RAI self-rule (civil)
me_r3d_m1g_t1_sc <- margins_summary(r3d_m1g_t1_sc, variables = "rai_autonomy", vcov = cluster.vcov(r3d_m1g_t1_sc, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 1)$cowcode)))
me_r3d_m1g_t1_sc$term <- "second-order maj."
me_r3d_m1g_t1_sf <- margins_summary(r3d_m1g_t1_sf, variables = "rai_autonomy", vcov = cluster.vcov(r3d_m1g_t1_sf, as.integer(subset(main_group, subset_analysis == 1 & int_grp_rel_dominant_g == 0)$cowcode)))
me_r3d_m1g_t1_sf$term <- "second-order min."
me_r3d_m1g <- rbind.fill(me_r3d_m1g_t1_sc, me_r3d_m1g_t1_sf)
me_r3d_m1g$term <- as.factor(me_r3d_m1g$term)
me_r3d_m1g$term = factor(me_r3d_m1g$term,levels(me_r3d_m1g$term)[c(2,1)])
me_r3d_m1g$model <- "RAI self-rule"
me_r3d_m1g$estimate <- me_r3d_m1g$AME
me_r3d_m1g$conf.low <- me_r3d_m1g$lower
me_r3d_m1g$conf.high <- me_r3d_m1g$upper
##f) RAI self-rule (communal)
me_r3d_m1d_t1_sb <- margins_summary(r3d_m1d_t1_sb, variables = "rai_autonomy", vcov = cluster.vcov(r3d_m1d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r3d_m1d_t1_sb$term <- "all"
me_r3d_m1d_t1_sb$model <- "RAI self-rule"
me_r3d_m2d_t1_sb <- margins_summary(r3d_m2d_t1_sb, variables = "rai_autonomy", at = list(included_excluded = c(0,1)), vcov = cluster.vcov(r3d_m2d_t1_sb, as.integer(subset(main_dyad,  subset_analysis == 1 & int_grp_rel_dominant == 1)$cowcode)))
me_r3d_m2d_t1_sb$term <- ifelse(me_r3d_m2d_t1_sb$included_excluded == 1, "included/excluded", "other")
me_r3d_m2d_t1_sb$model <- "RAI self-rule"
me_r3d_m13d <- rbind.fill(me_r3d_m1d_t1_sb, me_r3d_m2d_t1_sb)
me_r3d_m13d <- subset(me_r3d_m13d, term != "other")
me_r3d_m13d$term <- as.factor(me_r3d_m13d$term)
me_r3d_m13d$estimate <- me_r3d_m13d$AME
me_r3d_m13d$conf.low <- me_r3d_m13d$lower
me_r3d_m13d$conf.high <- me_r3d_m13d$upper
##g) combined (civil)
alt_aut_indices_plot_data_g <- rbind.fill(me_r3b_m1g,me_r3c_m1g, me_r3d_m1g)
alt_aut_indices_plot_data_g$model <- as.factor(alt_aut_indices_plot_data_g$model)
alt_aut_indices_plot_data_g$term= factor(alt_aut_indices_plot_data_g$term,levels(alt_aut_indices_plot_data_g$term)[c(2,1)])
alt_aut_indices_plot_data_g <- alt_aut_indices_plot_data_g[order(alt_aut_indices_plot_data_g$term),]
alt_aut_indices_plot_g <- dwplot(alt_aut_indices_plot_data_g,vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),dot_args = list(aes(colour = model,shape =model)), whisker_args = list(linetype = 1, aes(colour = model))) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("a) civil violence") +
  scale_color_manual(name = "Coefficient for:", values = c("#1b9e77","#d95f02","#7570b3"))+ scale_linetype_manual(name = "Coefficient for:", values = c(3,2,1))+ scale_shape_manual(name = "Coefficient for:", values = c(17,18,20))+ 
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of civil violence") + ylab("group type")+
  guides(shape = guide_legend(nrow=2,"model",reverse=TRUE), colour = guide_legend(nrow=2,"model",reverse=TRUE),linetype = guide_legend(nrow=2,"model",reverse=TRUE)) +
  scale_x_continuous(labels=scales::percent_format())
##h) combined (communal)
alt_aut_indices_plot_data_d <- rbind.fill(me_r3b_m13d,me_r3c_m13d,me_r3d_m13d)
alt_aut_indices_plot_data_d$model <- as.factor(alt_aut_indices_plot_data_d$model)
alt_aut_indices_plot_d <- dwplot(alt_aut_indices_plot_data_d,vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),dot_args = list(aes(colour = model,shape =model)), whisker_args = list(linetype = 1, aes(colour = model))) +
  theme_bw() + theme(plot.title = element_text(size=12)) + xlab("Coefficient estimate") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("b) communal violence maj./min. dyad") +
  scale_color_manual(name = "Coefficient for:", values = c("#1b9e77","#d95f02","#7570b3"))+ scale_linetype_manual(name = "Coefficient for:", values = c(3,2,1))+ scale_shape_manual(name = "Coefficient for:", values = c(17,18,20))+ 
  theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("change in predicted prob.\n of communal violence") + ylab("dyad type")+
  guides(shape = guide_legend(nrow=2,"model",reverse=TRUE), colour = guide_legend(nrow=2,"model",reverse=TRUE),linetype = guide_legend(nrow=2,"model",reverse=TRUE)) +
  scale_x_continuous(labels=scales::percent_format())
##i) combined
figurea18 <- grid_arrange_shared_legend(alt_aut_indices_plot_g, alt_aut_indices_plot_d, ncol=2, nrow=1)
ggsave(figurea18, file='../figures/figurea18.pdf', width = 14, height = 5.5, units="cm",dpi=1000)